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Abstract 

The identifiability of the two damping components of a Generalized Rayleigh Damping model is investigated through 
analysis of the continuum equilibrium equations as well as a simple spring-mass system. Generalized Rayleigh Damping 
provides a more diversified attenuation model than pure Viscoelasticity, with two parameters to describe attenuation 
effects and account for the complex damping behavior found in biological tissue. For heterogeneous Rayleigh Damped 
materials, there is no equivalent Viscoelastic system to describe the observed motions. For homogeneous systems, the 
inverse problem to determine the two Rayleigh Damping components is seen to be uniquely posed, in the sense that the 
inverse matrix for parameter identification is full rank, with certain conditions: when either multi-frequency data is available 
or when both shear and dilatational wave propagation is taken into account. For the multi-frequency case, the frequency 
dependency of the elastic parameters adds a level of complexity to the reconstruction problem that must be addressed for 
reasonable solutions. For the dilatational wave case, the accuracy of compressional wave measurement in fluid saturated 
soft tissues becomes an issue for qualitative parameter identification. These issues can be addressed with reasonable 
assumptions on the negligible damping levels of dilatational waves in soft tissue. In general, the parameters of a 
Generalized Rayleigh Damping model are identifiable for the elastography inverse problem, although with more complex 
conditions than the simpler Viscoelastic damping model. The value of this approach is the additional structural information 
provided by the Generalized Rayleigh Damping model, which can be linked to tissue composition as well as rheological 
interpretations. 
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Introduction 

The importance of damping models in elastography has become 
clearer in recent years as attenuation levels measured by 
elastographic imaging have been linked to diseases of the liver 
[1-6] and brain [7-13]. A number of methods have been 
proposed for reconstructing the Viscoelastic (VE) properties of 
soft tissue [14-18], including an iterative, nonlinear inversion 
method [19,20]. These methods have all targeted the development 
of images of the storage (G ) and loss (G ) modulus distributions 
within the tissue in question. Some have gone on to investigate the 
frequency dependent behavior of these two parameters [5,11], as 
well as multi-frequency reconstruction methods to improve the 
quality of the VE parameters reconstruction across a range of 
frequencies [18]. While these methods have already demonstrated 
the important role of tissue attenuation in differentiating tissue 
type and identifying lesions, linear VE provides a relatively 
simplified model for understanding the complex, non-linear 
attenuation observed in in-vivo tissue. 

Rayleigh Damping (RD), also known as proportional or 
Caughey damping, is a damping model with origins in numerical 
structural mechanics and is characterized by providing attenuation 
effects that are proportional to both elastic and inertial forces. As 
such, RD is a more diversified damping model than VE, where 
attenuation forces are related uniquely to elastic forces. From a 



numerical perspective, RD has the advantage that the damping 
matrix can be can be modally decomposed using the eigensystem 
developed from the undamped system, and has been shown to be 
useful in applications such as an absorbing boundary layer to 
remove spurious reflections in machine vibration and seismic 
models [21]. A rheological interpretation of RD can also be 
developed for weak to moderate damping levels [22]. 

A Generalized RD formulation has been developed previously 
for use in elastography imaging [23]. This time-harmonic 
formulation differs from the traditional RD configuration in the 
sense that damping effects are developed through complex valued 
shear modulus (//) and density (p) parameters, where the imaginary 
parts of ji and p are allowed to vary independently from their real 
counterparts. This is in contrast with the classical RD structural 
model, where the damping matrix is composed of scalar 
combinations of the mass and stiffness matrices. The difference 
is subtle, but important, as the use of the complex shear modulus 
value separates damping effects in distortional and dilatational 
waves, which will be seen to be critical for identifiability of the 
system. The use of Generalized RD in elastography is of interest 
mainly due to: a) the simplicity of the model, particularly in Finite 
Element (FE) formulations; b) the larger range of attenuation 
behavior the model is able to accommodate, where biological 
tissue is known to exhibit high levels of complex, non-linear 
damping; and c) the additional damping parameter provided by 
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the imaging process, which has shown sensitivity to material 
structure, such as the difference between gel and porous materials 
as well as cancerous and healthy breast tissue [24] . 

The goal of this paper is to investigate the conditions in which 
the parameters for a Generalized RD model can be identified in 
the elastography problem, based on measured motion data within 
the specimen. The reason to be concerned with identifiability in 
the RD case has to do with the simplified form of the classical RD 
reconstruction problem, which can be written as 

y,A+y 2 B=U, 

where complex valued A and B represent the stiffness and mass 
components of the structure. These parameters are to be identified 
based on measurements of complex valued U. While this case is 
only identifiable to a linear relation between A and B, 

A=iu-^B, 

Yi Yi 

it will be shown that this is in fact an over-simplification of the 
Generalized RD reconstruction problem that does not account for 
the dual modes of elastic wave transmission. Additionally, even in 
the single wave transmission case, it can be shown that the 
addition of multi-frequency information is sufficient to make the 
Generalized RD inverse problem uniquely posed. 

Analysis and Methods 

Generalized Time-Harmonic Rayleigh Damping 

Classical RD is a structural mechanics formulation based on the 
discretized dynamic model 



Mii + C« + Kw = 0, 



(1) 



with mass and stiffness matrices, M and K, and a damping matrix, 
C, of the form 



The elastography problem in a Generalized RD system is then 
to reconstruct the distribution of the elastic properties p R (x), Hj(x) 
and pi(x) from a measured displacement field, u(x). In general, 
for soft tissue elastography, the mass density, p R , can be assumed a 

kg 

priori as equal to that of water, i.e. 1000 — r. As a point of 

ifi 

comparison, the time-harmonic VE model provides damping 
effects due to a single elastic property, p { , such that the 
elastography problem in a purely VE system is to reconstruct 
p R (x) and p.j(x) from a measured displacement field, u(x). 

Defining the Generalized Rayleigh Damping 
Parameters. A physical interpretation of the Generalized RD 
parameters can be developed by considering the elastic equilib- 
rium of the time-harmonic system through Navier's equation, 
written as 



V-{fi{Vu + Vu T ))+V(XV-u) = -pw 2 u, 



(5) 



with complex valued fi and p defined as above in Eq. 3. In general, 
for soft tissue elastography, the system is considered nearly 
incompressible, with a real valued bulk modulus K ~0{\E9) and 



3 



(6) 



Dilatational attenuation can be considered, where 
K = K R + iKj, however the levels of damping in this compressional 
wave propagation will be negligible for higher frequencies where 
poroelastic effects are minimal [26]. 

Eq. 5 can be rewritten in the form of an undamped elastic 
operator YZ, a damping operator V and an inertial operator M, all 
acting on the time harmonic displacements, u. This gives 



JCu + iDu = —co 2 Mu, 



(7) 



with 



C = aM + /?K. 



(2) 



Km = V- (fi R (V« + Vu T ) ) + V(AjjV-m), 



(8) 



In contrast, the Generalized RD formulation considered in this 
work is a continuum interpretation of discretized classical RD, 
with damping effects proportional to shear and inertial forces. In 
the case of time-harmonic displacements, where 
u(x,i) = ^t{u(x)e' m '} for the complex valued u = u R + iui, the 
damping coefficients in Eq. 2 can be defined in terms of complex 
valued density, p = Pr + ipi, and shear modulus, p = p R -\- ipj, 
where [23] 

Pi = — ^ and n, = o)pp R . (3) 

We note here that the complex valued shear modulus 
components are directiy related to the storage and loss modulus, 
with p R = G and pj = G . The damping ratio, £, for an RD system 
is then given by [25], 




Vu = V fl u + V a u = [V- (fi, (VS + Vu T ) ) + V(//V-w)] + p^u (9) 
and 

Mu = p R u. (10) 

Substitution of Aj = wy). R and Eq. 3 into Eq. 9 and defining the 
time-harmonic velocity as v = ioM, gives the viscoelastic damping 
operator, 

C f v = iV ti u = V- (Pp R ( Vv + S/v T ) ) - V{yX R V- v) , (11) 

while substitution of Eq. 3 into Eq. 9 gives an inertial damping 
operator, 

C a y = iD a u = — ap R v. (12) 
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From Eqs. 11 & 12, the interpretation of the Generalized RD 
damping operator, C = Cfj + C a , is relatively straightforward. The 
component Cp represents traditional viscoelastic damping, where 
attenuation is related to strain rate. The component C x provides a 
second damping term directiy proportional to local velocity, as if 
the solid matrix were moving through a stationary viscous fluid. In 
the case of medical imaging, both of these components have a 
meaningful mechanical analogy to soft tissue structure. 

This Generalized RD formulation can be compared with the 
classical RD formulation by considering the case where y = fi, 
V-)3 = 0 and Vy = 0, so that 



Cpv=ifHCu, 



and 



(13) 



ICu + iVfjU = — (co 2 M + iT> a ,)u= — pco 2 u, 



(20) 



PR 



and then multiplying both sides by by — , which, after expansion 

of the derivative terms and combining p = Pr + ifij and 
a = Ar + iaj, gives 



jy w ^. wa+ (j^v 2 u+ (jyv^iwf 

+ ^^V(VS) r (21) 
+ (^f) (V; " )V '" + (^f ; ") W '" = ~ a2p R U - 



C a v= —W.M.U. 
We note that, in this case, Eq. 7 becomes 



-ia)Mu, 



(14) 



(15) 



which, as a single complex valued equation, cannot be solved for 
two independent parameters a. and /?. However, the conditions 
which lead to this case, notably the requirement that a and \i be 
essentially homogeneous, are not expected in in-vivo tissue, and Eq. 
6 is commonly used to define Xj, so that Eq. 9 becomes 



Cv = iVu = P 



V-(^(Vv+Vv r ))-V( 2 ^V-v 



-cr.p R v. 



(16) 



The implications of this particular case on the identifiability of 
the RD parameters will be seen in the following sections. 

Uniqueness of the Generalized Rayleigh Damping System 

A first step in analyzing the identifiability of the RD parameters, 
fi R , fij & pj, is to consider the existence of a VE system equivalent 
to the RD system but with a strictly real valued density of p R . If 
such a system exists, the RD parameters needed to define a 
particular motion field, u, axe not unique, and thus they cannot be 
identified without additional information. To develop the equiv- 
alence condition, we consider alternative operators for Eq. 7, so 
that we have 



K,u+iDu = K.u + iDu= — w 2 Mu, 



(17) 



where 



Making use of the product rule formulation 

a(db) = d(ab) — b(da), as well as effective moduli fi = ( — )p and 

P 

! = ( — )/. [24], Eq. 21 can be rewritten as 
P 



Ku+Vu- ( 



= — o) 2 Mu, 



■Vm- 



pV^iVuf 



P J (22) 



with the operators defined in Eqs. 18 & 19. We see from Eq. 22 
that Eq. 1 7 only holds for the case when the spatial derivative V 
is zero, i.e. when p is homogeneous or when p I = 0. In the case of 
heterogeneous density or pj^O, Eq. 22 indicates that some 
evidence of the intertial RD operator, D a , will be present in the 
motion field, u. We note that 



t,-t, j. iT. _( PrPr + PiPrPi \ i ,( PiP\-PrPrPi \ 
p-p R + ipj-\ 2 , + / 2 2 , 



V P 2 r + P 2 J V Pr + p] /' 



and 



l=l R+i l I= Cwr+^prP' ) +I f MpWA 

V Pr + Pi J V Pr + Pi J 



such that, even in the case where V^- 



0, the condition p I ¥ z 0 

ensures that fi R ^p R , fij^pj, Xr=£Xr and A/#A/. The presence 
of T> a in the system that generates u will thus effect measurements 
of both shear stiffness and viscosity made with a damping operator 

in the form of T>. 



JCu = Vip R (Vu + Vu T ))+v(i R V-uj, (18) 

and 

m=v\h{^u+Vu T ))+v{liS/-i^. (19) 



The form of K, and V can be determined by first moving D a to 
the right-hand-side of Eq. 7, to obtain 



Conclusions from Analysis of the Generalized System 

Eq. 22 shows that a direct VE equivalent to an RD system is 
Pr 

only possible when V — =0. This is an important result as, in 

general, the material property distributions observed within soft- 
tissue will be highly heterogeneous, and thus the RD and VE 
systems are, in general, not equivalent. The task still remains to 
correctiy identify the parameters of the RD system based on the 
measured motions. Eq. 15 has already indicated that this is 
generally possible so long as the tissue is heterogeneous and /? # y, 
i.e. dilatational wave attenuation is not at the same level as shear 
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wave attenuation. The exact conditions under which parameter 
identification is possible are analyzed in the following section. 

A Spring-Mass Analogy 

To explore the concept of parameter identification in Gener- 
alized RD elastography more closely, we start by considering a 
simple, locally homogeneous spring-mass system, where the spring 
stiffness, k, has a VE component, k , while the mass, m, has an 
inertial damping component, in . To focus on the case of 
elastography imaging, the spring system is constructed of three 
masses, connected by two springs, all in series, as shown in 
Figure 1. The reason for this arrangement is that typically, in 
elastography imaging, the applied forces, /, are unknown, and 
therefore cannot be called upon for the inverse problem of solving 
for the elastic parameters given the motions. Instead, only motion 
data is obtained, so in the two cases presented here we will 
consider the data as the axial displacements of the system, u\, U2 
and W3, measured at positions X\, X2 and X3. In addition, we will 
see that supplementary data can be obtained from measuring the 
corresponding transverse displacements, Vi, V2 and V3. Without 
forcing information, the trivial solution, m = k = m =k =0, is a 
possibility for the general inverse problem. To eliminate this 
solution, elastography methods typically assume the mass, m, is a 
known quantity. In soft-tissue imaging, this corresponds to the 
assumption that biological soft tissue has a density equal to that of 
water. 

A note on noise. In general, direct elastography inversion 
schemes such as those shown here are highly unstable, and the 
condition of the inversion matrix, Q, deteriorates rapidly with the 
presence of noise in the measured displacements. However, the 
purpose of this work is not to consider the impact of noise on the 
elastographic inverse problem. Any suitable approach to elasto- 
graphy imaging will have addressed the issue of the poor condition 
of the inverse solution matrix, and a number of well documented 
methods for resolving these issues have been presented in the 
literature. The choice of the direct inversion approach is made 
here because it allows explicit inspection of the inversion matrix to 
determine if the inverse problem is full rank. The sensitivity of 
these problems to noise in the data means that even if the inverse 
matrix is full rank, practical solutions from measured data requires 
filtering and regularization, despite the problem being "well 
posed". 



Case I: A 1D, Homogeneous System 

We start with the simplest of all cases, with a single mass, 
(m,m ), supported in series between two homogeneous springs, 
(k,k ), with steady displacements applied at both ends, (111,113), at 
radial frequency ft). 

Forward Problem. First, we consider the equation of motion 
for the interior mass at X2, U2 =f(m,m ,k,k ,u\,U2,co). This is 
obtained from considering the equilibrium equation for the 
system, given by 



dt 2 



d 2 U2 , / du2 , , , . , , 1 fdu 2 du\ 

m , „ +m — - +k(u2 — ui) — k(UT, — ii2) + k — — 

dt \dt dt 

1 / du 3 du 2 \ 
\dt dt J~ ' 



(23) 



which, with a time-harmonic assumption that the displacement 
will have the form u 2 = $l{A2e" 0 '}, (i = V — 1), leads to the 
equilibrium form 



— oj 2 m + icotn +k(u2 — ui) — k(u^ — U2) + io)k (u 2 — u\) 

— iwk (1/3 — u 2 ) = 0, 
with the solution 



(24) 



A 2 = 



(k + icok' j 



2k — a> 2 m + ico (m + 2k' 



TY^I+^S)- 



(25) 



Eq. 25 can be rewritten in the form 

A 2 = C(A 1 +A 3 ), 
where C has real and imaginary parts 



\k (2k - oj 2 m) + co 2 k' [m + 2fc') ] 
(2k — co 2 m) 2 + co 2 (m +2k') 2 



(26) 



and 



Q{C\- 



\o)k' (2k — a 2 m) ■ 



— wk(m +2k 



(2k — co 2 m) 2 + to 2 (m + 2k') 



(27) 



k, s 

•MAAd 




k,s 



VvW» 



Inverse Problem. The inverse problem in this case is ill 
posed, as can be deduced from the fact that the relationship 
between A2 and the driving conditions, Ai and A-}, is only governed 
by two numbers, 5J{C} and $i{C}. If we assume m is known, that 
leaves k, m and k to determine, with only two equations to use! 
The direct inverse problem can be written as 



Figure 1. A three mass spring-system as an elastography 
analogy. To eliminate the need for known forcing information in the 
elastography analysis, a three mass system is analyzed, where the 
displacements of masses m\ and m 3 are considered measured data and 
used to calculate system parameters. 
doi:1 0.1 371 /journal.pone.0093080.g001 



PLOS ONE I www.plosone.org 



4 



April 2014 I Volume 9 | Issue 4 | e93080 



Parameter ID in Rayleigh Damping Elastography 



?H{A 2 -A 1 -A 3 } -wQ{A 2 } -coQ{A 2 -Ai-A 3 } 
%{A 2 -A X -A 3 } oM{A 2 } oj^{A 2 -Ai-A 3 } 




o) 2 m&{A 2 } 
w 2 rr£${A 2 } 



The inverse solution matrix in this case, Q, is clearly only rank 
2, meaning we can compose the solution of any two parameters as 
functions of the remaining two. For example: 



k> = mQ{C} 2 w 2 + k' '<3{C}oj + m$t{C} 2 w 2 
23{C} 2 + 2K{C} 2 -K{C} 

and 



m (m,k ) = 

k'<3{C} 2 4+mwS{C}+k'yt{C} 2 4-k'$t{C}4 + k' 
23{C} 2 + 2K{C} 2 -5J{C} 



(30) 



Conclusions. The inverse problem in this case is not 
uniquely posed and the direct inversion matrix is rank deficient. 
There is no way to calculate a unique value of k and m given a 
single displacement measurement in a single dimension. We could 
consider adding an additional mass to the system to obtain more 
information, essentially adding another instance of Eq. 25, in the 
form 



(jc + iojk'^ 



Ai = - 2 V - (A 2 + A 4 ) 

2k — w L m + im[m +2k ) 



(31) 



Case II: A 1D System, with Multiple Frequencies 

The detection of the individual RD components requires 
additional, independent information in order to be uniquely 
posed. One option for obtaining this information is to consider 
measurements at different excitation frequencies, say co a and £0/,, 
such that 



■-$t{Ae"° a '} 



and 



u h = ft{Be' m b<}. 
The inverse problem system then becomes 



-ft) fl 3{^ 2 -^i-^3} 
mJR{A 2 -A\ -A 3 } 



1 

m 

(33) 



K{A 2 -Ai-A 3 } -co a Q{A 2 } 
%{A 2 -Ai-A 3 } WaR{A 2 } 
K{B 2 -B X -B 3 } -co b S{B 2 } —a>bSs{B 2 —Bi —B 3 } 
QiBi-B^-Bi} w,M{B 2 } w b $t{B 2 -B x -B 3 } 
{ w 2 a m&{A 2 } ) 
<j) 2 a m3{A 2 } 
0) 2 b mR{B 2 } 
, w 2 b mS{B 2 } 



We can examine the rank of this system by considering the 
determinant of the submatrix Q(l : 3, : ), whose numerator, for 
the case where A\ = A 3 = B\ = B 3 = j, is given by 



numerator{\Q(\ : 3, : )|} = 

OTft) a 3 ft) /) 2 (ft) a 2 — ft) A 2 ) (2k' +m'^j (k'm —kinj. 



(34) 



However, the additional measurement of A 3 provides no new 
information, as 



A 2 



(A l +A 3 ) (A 2 +A 4 )' 
and the expanded inversion matrix, Q, is still rank 2, with 



~${A 2 


-Ai 


-A,} - 


-wQ{A 2 } - 


w%{A 2 


-Ai 


-A 3 } 


<3{A 2 


-Ai 


-A,} 


oM{A 2 } w^{A 2 - 


Ai- 


-A 3 } 


$1{A 3 


-A 2 


~A 4 } - 


-coQ{A 3 } - 


coQ{A 3 


-A 2 


-A 4 \ 




-A 2 


-A 4 } 


0}^{A 3 } oj?H{A 3 - 


A 2 - 


-A 4 } 








' w 2 m^{A 2 } ' 














oj 2 mSs{A 2 } 


>. 












w 2 m>R{A 3 } 
_ w 2 m^s{A 3 } _ 










Eq. 34 has two real valued roots, (i> a = 0Jb and k =k — T (k and 

m 

m are positive reals, so the factor 2k +m will never be zero). The 
first of these roots reduces the system to the single mass, single 
frequency case described above. All additional information is 
identical to the original data, and the inversion matrix Q is once 
again rank 2. The importance of the second root is seen in 
considering the determinant of the submatrix Q(l; 2; 4, : ), which 
has the numerator 



numerator{\Q(\:2\4, : )|} = 

mo) a 3 o)i,(co a 2 — ft)/, 2 ) (2k—moj h 2 ) (jc m —lanj, 

where we see the common roots at m a = ft)/, and k 



(35) 



, m 
--k—. Both of 
m 

these conditions will thus reduce the inversion matrix to rank 2. 
The third root of 35, 2k — ma> b 2 , corresponds to undamped 
resonance at ft)/,, but due to the fact that this root does not appear 
in Eq. 34, this condition does not reduce the inverse problem 
matrix to rank 2. 

Conclusions. From Eqs. 34 & 35 we see that the multi- 
frequency inverse problem for RD systems is generally uniquely 
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posed, and, except in the case where k =k m r , Eq. 33 allows the 
determination of the two independent damping components k 
and m from data taken at two distinct frequencies. In practice, we 
note that material properties will often change with frequency, 
meaning that k, k , and even m should be expected to be different 
at frequencies w a and CO/,, thus making the inverse problem more 
complex than that posed in Eq. 33. This issue can be addressed by 
developing simple frequency dependency relationships for k, k , 
and even m , such as the power law relation, 8(a>) = 6(,of , and 
then adding additional frequency data, i.e. co c , to account for the 
additional parameters to be determined, at, a. k ' and ce m < . 



with the solutions 



and 



k + iwk 



2k — w 2 m + ico (m + 2k 



(41) 



( S + ICOS ) 

Bi=- 2 V / (B l+ B 3 ). 

2s — w L m + im[m +2s ) 



(42) 



Case III: A 2D System 

The frequency dependency of mechanical properties can be 
avoided in RD parameter reconstruction by considering a single 
frequency 2D wave propagation model for the spring-mass system 
shown in Figure 1. By allowing the mass to vibrate both axially, 
with displacement u = ^R{Ae"°'}, and transversely, with displace- 
ment v = $t{Be ,a "}, the model allows the propagation of both 
shear and longitudinal waves. From linear elastic theory, 
propagation of shear and longitudinal waves is governed by the 
elastic "moduli" s = fi and k = l + 2fi, respectively. For soft tissue 
imaging, it is common to make use of the fact that tissue is nearly 
incompressible due to its high water content. Thus, the bulk 
modulus, K, of the tissue can be specified at a very high numerical 
value (0(1 is 9)). The value of X can then be calculated from the 
relation 



X = K- 



2fi 

y 



which leads to the longitudinal wave "modulus" 



k=K + ^=K + 4 ^. 



(36) 



(37) 



To determine the VE component of k, it is assumed that K is 
real valued, so that k is determined by substituting the imaginary 
shear modulus component, Hj = Q{n}, into Eq. 37. We'll consider 
the case of complex valued K shordy. Thus, we have for the 
longitudinal wave attenuation 



= 0 + 



4hj 



4V 
3 



(38) 



Forward Problem. With the time-harmonic assumption 
noted above and the effective longitudinal wave modulus and 
attenuation given by Eqs. 37 & 38, we have the equilibrium 
conditions 



— oj 2 m + icom + k(u 2 — u\) — k(u 3 — u 2 ) J \-iwk (u 2 —u\) 



-iwk («3 — u 2 ) = 0, 



(39) 



and 



Inverse Problem. The inverse problem system then becomes 



?H{A 2 -A l -A i } -co a S{A 2 } -4fL${A 2 -Ai-A3} 



4<»a. 



m 
(43) 



-c^a 2 -Ai-A 3 } wJSt{A 2 } -±fSt{A 2 -Ai-A 3 } 

$t{B 2 -Bi-B 3 } -co h Z{B 2 } -(O b S{B 2 -Bi-B 3 } 
Q{B 2 -By-Bi} oj h K{B 2 } w h ^{B 2 -B l -B,} 
' co 2 a m^{A 2 } - m{A 2 -Ai-A 3 }' 
co 2 a mQ{A 2 }-K^{A 2 -A l -A ? ,} 
oj 2 h mR{B 2 } 
w 2 h mSs{B 2 } 



We can examine the rank of this system by considering the 
determinant of the submatrix Q(l : 3, : ), whose numerator, for 

the case where A\ =A 3 = B\ =B 3 = -, is given by 



numerator{\Q(\ : 3, : )|} = 
12^w 5 (m 2 +m 2 co 1 



[m 1 +m 2 co 2 ) (m'+2s^ 



(44) 



Eq. 44 has no real roots, given that s >0 and in >0 for RD 
systems. We note also that the determinant of the submatrix 
Q(l; 2; 4, : ) has the numerator 



numerator{ \ Q(l; 2; 4, : )|} = 



\2Ka> 



1 



{nl 1 +m 2 oj 2 ^ (2s- 



(45) 



with a root at i= -mco , which corresponds to the undamped 

resonance case in the shear displacements, v. However, as this root 
only appears in Eq. 45, this condition does not reduce the inverse 
problem matrix to rank 2. 

Conclusions. From Eqs. 44 & 45 we see that the 2D problem 
for RD systems is generally uniquely posed, and Eq. 43 allows the 
determination of the two independent damping components k 
and m from shear and longitudinal wave data. 



— co 2 m + icom + s(v 2 — v\ ) — .s(v 3 — v 2 ) + icos (v 2 — vi ) 
-iws (vi - v 2 ) =0, 



(40) 
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Case IV: A 2D System, with Relative Damping 
Components 

The 2D problem above can be generalized for a system where 
attenuation occurs in both dilatational and distortional wave 
propagation. In this case, the shear and longitudinal stiffnesses, s 
and k, as well as their VE components, S and k , can be direcdy 
related by factors a r and a,-, such that 



whose numerators, for the case where A\ = A 3 = B\ = B 3 = -, are 
given by 

numerator{\Q(\ : 3, : )|} = 

, (2<x r m s 1 — 2<Xi>n' co 2 s 2 — atm 2m 2 s — a,'m 2 co 4 /\ (51) 

-mco («,-« r ) 

\ + 2a,mtu z i s + 2a r mw^s s 



and 



and 



k =tXjS . 



numerator{\Q(\;2;4, : )|} = 



-ma> 6 (otj — a r ) 



a r »! j — 2i r ms +2otitn s s + 2ot r m s s 
+ 2a i mcoi 2 s' 2 + a r m 2 w 2 s 



(52) 



Forward Problem. With the time-harmonic assumption, we 
have the equilibrium conditions 



-co 2 m + ia>m + oc r s(u 2 —ui) — cc r s(u 3 — u 2 ) + wjaiS (ii2 — u\) 



-itoa.jS (u 3 — u 2 ) =0, 



and 



(46) 



respectively. We note that both numerators share the factor a; — oc r , 
meaning that the inverse problem becomes ill posed in the case 

k' k 
where — = — . 
s s 

Conclusions. From Eqs. 5 1 & 52 we see that the 2D problem 
for RD systems is generally uniquely posed, except for the case 
where the damping ratios for distortional and dilatational wave 
attenuation are equal, i.e. whenever 



— co 2 m + icom + s(v 2 — vi) — .s(v 3 — v 2 ) + icos ( v 2 — vi ) 

-KOS (V3-V2)=0, 

with the solutions 



(47) 



k_ 
k 



s 
s 



(53) 



Equivalent to Eq. 4, the damping ratio, £, for wave propagation 
in a mode governed by stiffness T is given by the relation 



CC r S+lOJ0LiS 



2rx r s — 03 2 m + w} (m + 2a 



tt(Ai+A 3 ) (48) 



and 



2s — co 2 m+ico(m'+ 2s ) 



_ (49) 



Inverse Problem. The inverse problem system then becomes 



a. r ^.{A 2 -A x -A 3 } -oj a Q{A 2 } -w a a i ^s{A 2 - A A - A 3 } 
4 

a r -^{A 2 -A\-A 3 } coJR{A 2 } co a a.^{A 2 -Ai-A 3 } 

^{B 2 -B x -B 3 } -w b Q{B 2 } -a b SS{B2-Bi-B3} 

Q{B 2 _ Bl -B 3 } w h Vt{B 2 } w b %t{B 2 -B x -Bi} 



co 2 a m$t{A 2 } * 
(i> 2 a mSs{A 2 } 
o) 2 h m&{B 2 } 
k a) 2 mQ{B 2 } 



(50) 



r 1 / T m I 

2 \ T w 



where x is the VE attenuation in the given propagation mode {i.e. 
either s for shear waves or k for longitudinal waves). In this case, 
we see that the condition for which the inverse problem becomes 
ill posed, given by Eq. 53, is equivalent to the damping ratio for 
wave attenuation for the shear and longitudinal waves being equal. 
In this case, no new information is added to our inverse problem 
system by considering the two different wave propagation cases, 
and the problem becomes rank deficient. 

Case V: A 2D System, with Dilatational Damping 

The generalized case discussed above breaks down in the case 
where a r = a,-. In the case of damped dilatational waves, the 
coefficients <i r & a,- are defined as 



4s (K 4\ 



and 



k =a/S =K + 



As 



K 4^ 



We can examine the rank of this system by considering the 
determinant of the submatrices Q(l : 3, : ) and Q(l;2;4, :), 



The condition cc r = a, is then equivalent to 
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AT 
A 7 



s 

s ' 



i.e. the condition of equivalent damping ratios for distortional and 
dilatational wave propagation. This singularity condition can be 

eliminated by considering K as a separate unknown in the RD 
inverse problem (we can maintain the "nearly incompressible" 
condition by assuming K ~0(\E9)). 

Forward Problem. With the time-harmonic assumption, we 
have the equilibrium conditions 



4s 
1 



4?Y 



co m + iiom + ( K+ — ] («2 — Mi)— ( K+ —J (M3 — 112) 
4s \ , , / / 4s 



+ ioj\K + — Uu 2 -ui)-uolK +— |(m 3 -m 2 )=0, 



(54) 



We can examine the rank of this system by considering the 
determinant of the full matrix Q, whose numerator, for the case 

where A\ = At, = B\ = B^ = —, is given by 



numerator {\Q\} = 

9a) 6 (m' 2 +m 2 w 2 ) ( 36K " ma? + 36K~ Km 0/ + 96K'm afs 



- 96K 2 mw 4 s + 36K' K 2 mco 2 -9K m' 2 mw 4 +4%Km2w 2 s 

- U8K in oj 2 s s-9K m^co 6 + 48K m 2 o) 4 s + 64K mco 4 / (59) 

- 64K mco 2 s 2 + 36K 3 m + 96K 2 m s + 96K 2 mw 2 s 

- 9Km l co 2 -4%Km 2 (a 2 s - 9Kmm 2 w 4 - 64Km co 2 s' 2 



+ (AKm s 2 -4%Km 2 m 4 s + W&Kmco 2 s s ) , 



and 



-a m + team +s(\>2 — Vi) — s(y$ — V2) + /co.s (v'2 — Vi) 

-IWS (t'3 — V2) =0, 



with the solutions 



(55) 



whose roots, given as values of AT =f\s, s , m J , are complex 

enough to essentially eliminate the chance of them occurring for 

reasonable values of 5, s & m in soft tissue. 

Conclusions. Eq. 59 indicates that the general RD inverse 
problem, even in cases with dilatational wave attenuation, is 
generally uniquely posed once we include K among the 
parameters for reconstruction. While an interesting theoretical 
result, the practicality of obtaining meaningful results from 
dilatation measurements in Elastography is severely limited, given 
their long wavelengths and the the susceptibility of these 
measurements to noise. 



K +J ) +l co(K + - 



2( K+j\ -w 2 m + ko(m' +2( K' + ^ 



7YV (^i+^ 3 )(56) 



and 



Bi- 



2s — oi 2 m + 10) (in + 2s ) 



TV {B1+B3) 



(57) 



Inverse Problem. The inverse problem system then becomes 



4 4a)„ 

^U{A 2 -A, -A,} -w a ^{A 2 } --^H{A2-Ai -A,} -^{A 1 -Ai -A 3 
4 4<u, 

-HAl-Ai-A,} coMAl] -^mAl-Ai-Ai} K{A 2 -Ai-A,} 



K{B 2 -B I -B 3 } -w b <s{B 2 } -m b Q{B 2 -Bi-B,} 
3{B 2 -B,-B } } co s »{S 2 } <oiM{B 2 -B,-B 3 } 
( a$m'gt{A 2 }-K , Sl{A 2 -A, -A,} ' 
o? a m^{A 2 ] -KS{A 2 -A, -A,} 
a>lm%t{B 2 } 



(58) 



Conclusions from Spring-Mass Analysis 

A summary of the conclusions from the above cases is given 
below: 

• From Case I we see that the simple problem of deducing from 
k, k and m from the ID displacement of a single mass at a 
single frequency is impossible. This is not a surprise, as the data 
here consists of a single complex number, i.e. two measure- 
ments cannot produce three parameters. The inverse problem 
matrix, Q, is rank 2. 

• From Case II we see that, in general, the addition of data taken 
at another frequency makes the RD inversion problem 

uniquely posed. The inverse problem matrix is rank 3 except 
/ m 

in the special case where k =k — r , which corresponds to a 

m 



damping ratio t, ■ 



m 2 + (mf 



From Case III, and more generally Case IV, we see that the 
addition of shear wave data from the transverse direction 
makes the problem uniquely posed, so long as the attenuation 
differs between the shear and longitudinal wave propagation. 

From Case V, we see that the general problem of RD 
parameter reconstruction becomes uniquely posed if we 
include the dilatational attenuation, K , as an unknown 

k s 

parameter, even in the case where — = — 
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We see that the spring-mass RD reconstruction problem is 
generally uniquely posed with additional measurements, either 
from the propagation of waves along a different mode, or from 
additional frequency information. In practice, these two informa- 
tion sources contain their own challenges. In the case of 
measurements at additional frequencies, the frequency dependen- 
cy of the parameters comes into play, and in itself requires 
additional information (as well as good models for this frequency 
dependency) in order to pose a reasonable reconstruction problem. 
In the case of deducing information from longitudinal waves in 
elastography, the problem arises from the nearly incompressible 
nature of soft tissue, where the longitudinal wavelength becomes 
very long (D(lm) in comparison to the size of the objects being 
measured O(10cot). Accurately characterizing these wavelengths 
within the medium is thus highly susceptible to noise. 

Discussion 

Overall, we can say that the RD parameters are identifiable, 
given certain conditions and assumptions on the damping 
behavior of the elastic material. To start with, in a region of 
heterogenous material properties, where the term V y is unlikely 
to disappear, analysis of the generalized RD system shows that it 
has no purely VE equivalent, and thus there is a valid reason to 
consider RD reconstruction. Next, to determine the RD param- 
eters, the assumptions required for identification are not partic- 
ularly onerous, specifically, the idea that the attenuation of 
dilatational waves is of a different order to the attenuation of shear 

References 

1. Rouviere O, Yin M, Dresner M, Rossman P, Burgart L, ct al. (2006) MR 
elastography of the liver: Preliminary results. Radiology 240: 440^448. 

2. Klatt D, Asbach P, Rump J, Papazoglou S, Somasundaram R, et al. (2006) In 
vivo determination of hepatie stillness using steady-state free preeession magnetie 
resonance elastography. Investigative Radiology 41: 841—848. 

3. Huwart L, Sempoux C, Salameh N, Jamart J, Annet L, et al. (2007) Liver 
fibrosis: Noninvasive assessment with mr elastography versus aspartate 
aminotransferaseto-platelet ratio index. Radiology 245: 458-466. 

4. Klatt D, Hamhaber U, Asbach P, Braun J, Sack I (2008) Noninvasive assessment 
of the rheologieal behavior of human organs using multifrequcncy mr 
elastography: a study of brain and liver viscoelasticity. Physics in Medicine 
and Biology 52: 7281-7294. 

5. Klatt D, Frcidrich C, Korth Y, Vogt R, Braun J, et al. (2010) Viscoelastic 
properties of liver measured by oscillatory rheomctry and multifrequcncy 
magnetic resonance elastography. Biorheology 47: 133—141. 

6. Venkatcsh S, Yin M, Ehman R (2013; Magnetic resonance elastography of liver: 
technique, analysis, and clinical applications. Journal of Magnetic Resonance 
Imaging 37: 544-555. 

7. Vappou J, Willinger R (2006) Dynamic viscoelastic shear properties of soft 
matter by magnetic resonance elastography using a low-field dedicated system. 
Journal of Rhcology 50: 531-542. 

8. Sack I, Bcierbach B, Hamhaber U, Klatt D, Braun J (2007) Non-invasive 
measurement of brain viscoelasticity using magnetic resonance elastography. 
NMR in Biomedicine 21: 265-271. ' 

9. Green M, Bilston L, Sinkus R (2008) In vivo brain viscoelastic properties 
measured by magnetic resonance elastography. NMR in Biomedicine 21: 755- 
764. 

10. Wucrfcl J, Paul F, Beierbach B, Hamhaber U, Klatt D, et al. (2010) Mr- 
elastography reveals degradation of tissue integrity in multiple sclerosis. 
Ncuroimagc 49: 2520-2525. 

11. Clayton E, Garbow J, Bayly P (2011) Frequency-dependent viscoelastic 
parameters of mouse brain tissue estimated by mr elastography. Physics in 
Medicine and Biology 56: 2391-2406. 

12. Zhang J, Green MA, Sinkus R, Bilston LE (2011) Viscoelastic properties of 
human cerebellum using magnetic resonance elastography. Journal of 
Biomechanics 44: 1909-1913. 

13. Murphy M, Curran G, Glascr K, Rossman P, Huston J III, ct al. (2012) 
Magnetic resonance elastography of the brain in a mouse model of alzheimcr's 
disease: initial results. Magnetic Resonance Imaging 30: 535-539. 



waves is quite reasonable and easy to justify. The real 
issue, however, is the dependence of the reconstruction on 
measurements of the dilatational wave component itself, which 
are easily corrupted by the presence of noise. One way to alleviate 
this problem is by introducing multiple frequency data into the 
reconstruction problem, which, except in special circumstances, 
renders the RD identification problem uniquely posed. Multi- 
frequency reconstruction introduces its own set of complications 
however, due to the frequency dependence of the parameters in 
the elastic equilibrium equations. In short, RD parameter 
identification is possible, but not a simple affair. The intrigued 
reader asking "why bother?" is encouraged to consult the results 
presented in [24] for some evidence of the potential value of RD 
parameter data in elastography imaging. 

Acknowledgments 

The author would like to thank Ralph Sinkus for asking the questions that 
eventually led to this work, and Paul Barbone for his insight in formulating 
the response. 

Elijah E.W. Van Houten is a member of the FRSQjfunded Centre de 
recherche clinique Etienne-Le Bel. 

Author Contributions 

Conceived and designed the experiments: EVH. Performed the experi- 
ments: EVH. Analyzed the data: EVH. Contributed reagents/materials/ 
analysis tools: EVH. Wrote the paper: EVH. 



14. Sinkus R, Sicgmann K, Tanter M, Xydcas T, Catheline S, et al. (2005) 
Viscoelastic shear properties of in vivo breast lesions measured by MR 
elastography. Magnetic Resonance Imaging 23: 159—165. 

15. Eskandari H, Salcudean S, Rohling R (2008) Viscoelastic parameter estimation 
based on spectral analysis. IEEE Trans Ultrason Ferroelcctr Frcq Control 55: 
1611-1625. 

16. Bayly P, Clayton E, Genin G (2012) Quantitative imaging methods for the 
development and validation of brain biomechanics models. Annual Review of 
Biomedical Engineering 14: 369-396. 

17. Lcclcrc G, Dcbcrnard L, Foucart F, Robert L, Pcllcticr K, ct al. (2012) 
Characterization of a hyper-viscoclastic phantom mimicking biological soft tissue 
using an abdominal pneumatic driver with magnetic resonance elastography. 
Journal of Biomechanics 45: 252-257. 

18. Papazoglou S, Hirsch S, Braun J, Sack I (2012) Multifrequcncy inversion in 
magnetic resonance elastography. Physics in Medicine and Biology 57: 2329- 
2346. 

19. Doylcy M, Perrcard I, Patterson A, Weaver J, Paulsen K (2010) The 
performance of steady-state harmonic magnetic resonance elastography when 
applied to viscoelastic materials. Medical Physics 37: 3970-3979. 

20. Mcgarry M, Van Houten E, Johnson C, Georgiadis J, Sutton B, et al. (2012) 
Multiresolution mr elastography using nonlinear inversion. Medical Physics 39: 
6388-6396. 

21. SemblatJ, Lend L, Gandomzadch A (201 1) A simple multi-directional absorbing 
layer method to simulate elastic wave propagation in unbounded domains. 
International Journal For Numerical Methods In Engineering 85: 1543-1563. 

22. Scmblat J (1997) Rheologieal interpretation of rayleigh damping. Journal of 
Sound and Vibration 206: 741-744. 

23. McGarry M, Van Houten E (2008!; Use of a rayleigh damping model in 
claslography. Medical and Biological Engineering and Computing 46: 759-766. 

24. Van Houten E, Viviers D, McGarry M, Perrifiez P, Perreard I, et al. (2011) 
Subzonc based magnetic resonance elastography using a rayleigh damped 
material model. Medical Physics 38: 1993-2004. 

25. Cook R, Malkus D, Plcsha M, Witt R (2002) Concepts and Applications of Finite 
Element Analysis. John Wiley & Sons, Inc., 4th edition. 

26. Mcgarry M, Pattison A, Van Houten E, Johnson C, Sutton B, ct al. (2012) 
Analysis of viscoelastic and poroelastic behavior in mre. In: Proc. of ISMRM 
2012. ISMRM, p. 2520. 



PLOS ONE | www.plosone.org 



9 



April 2014 | Volume 9 | Issue 4 | e93080 



